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Abstract We propose a multiscale mechanobiological 
model of bone remodelling to investigate the site-specific 
evolution of bone volume fraction across the midshaft 
of a femur. The model includes hormonal regulation 
and biochemical coupling of bone cell populations, the 
influence of the microstructure on bone turnover rate, 
and mechanical adaptation of the tissue. Both micro¬ 
scopic and tissue-scale stress/strain states of the tissue 
are calculated from macroscopic loads by a combination 
of beam theory and micromechanical homogenisation. 

This model is applied to simulate the spatio-temporal 
evolution of a human midshaft femur scan subjected to 
two deregulating circumstances: (i) osteoporosis and 
(ii) mechanical disuse. Both simulated deregulations led 
to endocortical bone loss, cortical wall thinning and 
expansion of the medullary cavity, in accordance with 
experimental findings. Our model suggests that these 
observations are attributable to a large extent to the 
influence of the microstructure on bone turnover rate. 
Mechanical adaptation is found to help preserve intra- 
cortical bone matrix near the periosteum. Moreover, 
it leads to non-uniform cortical wall thickness due to 
the asymmetry of macroscopic loads introduced by the 
bending moment. The effect of mechanical adaptation 
near the endosteum can be greatly affected by whether 
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the mechanical stimulus includes stress concentration 
effects or not. 
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1 Introduction 

Bone is a biomaterial with a complex hierarchical struc¬ 
ture characterised by at least three distinct length scales: 
(i) the cellular scale (10-20 /xm); (ii) the tissue scale (2- 
5 mm) and (hi) the whole organ scale (4-45 cm) [Rho 
et al (1998); Weiner and Wagner (1998)]. Several in¬ 
teractions exist between these scales, which affect bone 
remodelling, bone material properties and bone struc¬ 
tural integrity. The activity of bone-resorbing and bone¬ 
forming cells during bone remodelling leads to changes 
in material properties at the tissue scale which subse¬ 
quently affect the distribution of loads at the structural, 
whole organ scale (Figure 1). Besides, changes in bone 
shape and microarchitecture modify the stress/strain 
distribution and bone surface availability, which provide 
mechanical and geometrical feedbacks onto the bone 
cells and, eventually, affect bone remodelling [Martin 
(1972); Lanyon et al (1982); Frost (1987)]. Due to the 
complexity of these interactions, the interpretation of 
experimental data at a single scale is difficult. Predict¬ 
ing the evolution of multifactorial bone disorders, such 
as osteoporosis, necessitates a comprehensive modelling 
approach in which these multiscale interactions are con¬ 
sistently integrated. 

Various mathematical models of bone remodelling 
have been proposed in the literature. Biomechanics mod¬ 
els estimating tissue-scale stress and strain distribu¬ 
tion from musculoskeletal models and average material 
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(a) Couplings in the bone 
remodelling process 


(b) Femur (c) Midshaft cross-section Representative volume of 

the tissue at r 


Biochemical couplings 



Osteoblasts Osteoclasts 



£tissue(^^ ?) 



Fig. 1 Multiscale representation of bone, (a) Scheme of the couplings in the bone remodelling process; (b) Femur bone geometry 
(organ scale); (c) Midshaft cross section depicting coordinate axes and the sectional forces used for beam theory (tissue to organ 
scales); (d) Representative volume element (RVE) of cortical bone used to define bone cell densities, bone volume fraction, and 
specific surface (cellular to tissue scales). 


properties, such as bone density, are often used in con¬ 
junction with remodelling algorithms based on Wolff’s 
law. These remodelling algorithms locally increase or 
decrease bone density depending on the tissue’s mechan¬ 
ical state [Carter and Hayes (1977); Carter and Beaupre 
(2001); Fyhrie and Carter (1986); Weinans et al (1992); 
Van der Meulen et al (1993); Pettermann et al (1997)]. 
Such models may also include damage accumulation 
due to fatigue loading and damage repair [Prendergast 
and Taylor (1994); McNamara and Prendergast (2007); 
Garcia-Aznar et al (2005)]. Other models focus at the 
microstructural scale (/rm to mm) and describe the evo¬ 
lution of trabecular bone microarchitecture through re¬ 
sorption and formation events at the bone surface in¬ 
duced by the local mechanical state [Huiskes et al (2000); 
Ruimerman et al (2005); Van Oers et al (2008); Christen 
et al (2012, 2013)]. Most of these mathematical models 
focus on the biomechanical aspects of bone remodelling 
and do not consider hormonal regulation or biochemical 
coupling between bone cells. 

In this paper, we propose a novel multiscale mod¬ 
elling approach of bone remodelling combining and ex¬ 
tending several mathematical models into a consistent 
framework. This framework enables (i) the considera¬ 
tion of biochemical and cellular interactions in bone 
remodelling at the cellular scale [Lemaire et al (2004); 
Pivonka et al (2008); Buenzli et al (2012); Pivonka et al 
(2013)], (ii) the evolution of material properties at the 
tissue scale based on bone cell remodelling activities reg¬ 
ulated by mechanical feedback [Scheiner et al (2013)] 
and bone surface availability [Pivonka et al (2013); Buen¬ 
zli et al (2013)], and (iii) the determination of the stress/ 
strain distributions from the tissue scale to the mi¬ 


crostructural scale by a combination of generalised beam 
theory and micromechanical homogenisation [Hellmich 
et al (2008); Scheiner et al (2013); Buenzli et al (2013)]. 


This modelling approach is applied to simulate the 
temporal evolution of a human femoral bone at the 
midshaft (Figure 1), subjected to various deregulating 
circumstances such as osteoporosis and changes in me¬ 
chanical loading. An initial state of normal bone remod¬ 
elling is first assumed, in which the tissue across the 
midshaft cross section remodels at site-specific turnover 
rates without changing its average material properties. 
Osteoporosis is then simulated by hormonal changes 
deregulating the biochemical coupling between osteo¬ 
clasts and osteoblasts. These hormonal changes are cal¬ 
ibrated so as to reproduce realistic rates of osteoporotic 
bone loss. The strength of the resorptive and forma¬ 
tive responses of bone cells to mechanical feedbacks are 
calibrated so as to reproduce rates of bone loss and re¬ 
covery in cosmonauts undertaking long-duration space 
flight missions. A scan of a femur cross section is used 
as initial condition for our simulations. This illustrates 
the potential of our modelling approach to be used as 
a predictive, patient-specific diagnostic tool for estimat¬ 
ing the deterioration of bone tissues. Here, we use the 
model to investigate the interplay between geometrical 
and mechanical feedbacks in inducing site-specihc bone 
loss in osteoporosis, which is characterised by endocorti- 
cal bone loss, cortical wall thinning, and the expansion 
of the marrow cavity [Feik et al (1997); Bousson et al 
(2001); Zebaze et al (2010)]. 
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2 Description of the model 

Figures 1 and 2 summarise the general approach of our 
model. We consider a portion of human femur near 
the midshaft. This portion of bone is assumed to carry 
loads corresponding to a total normal force N and to¬ 
tal bending moment M (Figure 1(c)). These loads are 
distributed unevenly across the midshaft cross section 
depending on the site-specific bone microstructure, par¬ 
ticularly on the cortical porosity [Zebaze et al (2010); 
Buenzli et al (2013)]. This load distribution determines 
a site-specific mechanical stimulus which is sensed and 
transduced by bone cells (Figure 2(a)). This mechanical 
feedback is incorporated in a cell population model as 
biochemical signals leading to changes in the balance 
between osteoclasts and osteoblasts (Figure 2(b)). In 
addition, microstructural parameters such as bone vol¬ 
ume fraction (/bm) and bone specific surface influence 
the propensity of bone cells to differentiate and become 
active [Martin (1984); Lerebours et al (2015)]. This ge¬ 
ometrical feedback is included in the cell population 
model via a dependence of the bone turnover rate on 
the bone volume fraction. The activities of osteoclasts 
and osteoblasts modify the tissue microstructural pa¬ 
rameters (bone volume fraction, bone specific surface), 
which in turn induces changes in the load distribution 
(Figure 2(c)). In the following, we introduce in more 
detail the multiple scales and related variables involved 
in this model workflow. Table 1, in Appendix A.2, lists 
all the parameters of the model. 

2.1 Load distribution from the organ scale to the 
cellular scale 

Loading is composed of body weight and muscle forces 
exerted onto bone via tendons and direct action of mus¬ 
cles. These forces can be calculated from bone shape, 
muscle and tendon attachment, and gait analysis data 
using musculoskeletal models [Lloyd and Besier (2003); 
Viceconti et al (2006); Martelli et al (2014)]. Continuum 
mechanics provides the link between external forces ex¬ 
erted onto a structure, and the strain and stress distri¬ 
bution in the structure [Salencon (2001)]. 

Tissue-scale properties within the framework of con¬ 
tinuum mechanics are average mechanical properties 
over microstructural material phases and pores (pre¬ 
sented in detail in the next sections). The corresponding 
tissue-scale stresses and strains may significantly deviate 
from the microscopic, cellular scale, stresses and strains 
acting in the different material phases composing the 
tissue due to so-called strain and stress concentration ef¬ 
fects [Zaoui (2002); Hill (1963)]. Microscopic stress and 
strain distributions in the bone matrix are likely to be 


sensed directly by bone cells, particularly by osteocytes 
[Scheiner et al (2013)]. However, as osteocytes form an 
extensive interconnected network [Marotti (2000); Buen¬ 
zli and Sims (2015)], they may also sense larger scale 
stress and strain distributions. We will let either the 
tissue-scale or the microscopic mechanical state of bone 
act onto the bone cells to investigate how this influences 
the site-specific evolution of bone tissue microstructures. 

In the following, we present first how stress and 
strain distributions can be calculated at the tissue scale 
using beam theory. We then present how these tissue- 
scale stress and strain distributions are employed as 
site-specific loading boundary conditions to the contin¬ 
uum micromechanical model of Hellmich et al (2008) for 
the calculation of microscopic stress and strain distribu¬ 
tions effective at the cellular level. 

Determination of tissue-scale stress and strain distribu¬ 
tions based on beam theory 

The continuum mechanical field equations allow the cal¬ 
culation of tissue-scale strain and stress distributions in 
bone. Given that the length of the femur L is signifi¬ 
cantly larger (45-50 cm) than its diameter D (3-5 cm) 
at the midshaft (Figure 1(b)) the continuum mechanical 
field equations can be approximated using beam theory 
formulated for small strains and generalised to materials 
of non-uniform properties, an approach we have used 
previously in Buenzli et al (2013). 

This approach requires the knowledge of the total 
external forces, i.e., the normal force N and the bend¬ 
ing moment M carried by the femur cross section. TV 
and M can be estimated for different physical activ¬ 
ities by using musculoskeletal models [Vaughan et al 
(1992); Forner-Cordero et al (2006)]. In our simulations 
we take constant values for TV and M comparable with 
the maximum ground reaction force and knee and hip 
moments that occur during a gait analysis, estimated 
as; TV = {N^, 0,0), = -700 N, and M = TH m, M 

= 50 Nm, where rh, is a unit vector along the antero¬ 
posterior axis of the cross section determined from the 
micro-radiograph [Vaughan et al (1992); Forner-Cordero 
et al (2006); Ruff (2000); Gordey and Gautier (1999)] 
(Figure 5(c)). The x-axis is the femur’s longitudinal axis 
and {y,z) is the plane transverse to x at the midshaft 
(Figure l(b)-(c)). 

Tissue-scale mechanical properties correspond to spa¬ 
tial averages over a so-called representative volume el¬ 
ement (RVE) of the tissue. In cortical bone, an appro¬ 
priate tissue RVE is of the order of 10 x 2 x 2 mm^, a 
size large enough to contain a large number of pores, 
but small enough to retain site-specific information and 
to not be influenced by macroscopic features such as 
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Organ scale Tissue and cellular scales 



Fig. 2 Flow chart of bone remodelling simulations taking into account (a) the global mechanical loading, (b) the bone cell 
population model, and (c) the bone material and geometry adaptation. 


overall bone shape [Hill (1963); Zaoui (2002)]. We de¬ 
note by t) the local bone tissue stiffness tensor 

defined at the RVE scale, where r denotes the location 
in bone of the RVE (Eigure 1(c)) and the dependence 
on time t reflects the fact that bone remodelling may 
modify the local mechanical properties of the tissue. 
This tissue-scale stiffness tensor is assumed to relate 
the tissue-scale stress tensor and strain tensor 

^tissue pointwise according to Hooke’s law: 


’(r,t) =C‘"®“(r,t) :e*‘"™®(r,t). 


( 1 ) 


Beam theory is based on the so-called Euler-Bernoulli 
kinematic hypothesis, which asserts that material cross 
sections initially normal to the beam’s neutral axis re¬ 
main planar, undeformed in their own plane, and normal 
to the neutral axis in the beam’s deformed state [Timo¬ 
shenko and Goodier (1951); Bauchau and Craig (2009); 
Hjelmstad (2005)]. These assumptions are expected to 
be well satisfied near the femoral midshaft under small 
deformations generated by bending and compression 
or tension. Furthermore, no shear force, torsional loads 
or twisting along the beam axis are assumed. These 
assumptions, Eq. (1), and the fact that bone is an or¬ 
thotropic material [Hellmich et al (2004)] imply that 
the only nonzero components of the stress tensor are 
the normal stresses "®, = 

and cr“f’^® = where the nor¬ 
mal stresses cr**®®"® and are induced by compres¬ 

sion or tension along the beam axis x by the Poisson 
effect^ (see Buenzli et al (2013) for more details). The 


^ The stress components and do not partic¬ 

ipate directly to the transfer of the resultant force N and 
resultant bending moment M across the bone cross section, 
however, they are accounted for in the calculation of the tissue- 
scale strain energy density 


Euler-Bernoulli hypothesis implies that the tissue strain 
tensor reduces to the single non-zero component 
and that: 

= eiit) - K3it)y + K2it)z, ( 2 ) 

where Si is the sectional axial strain, and K 2 and K 3 are 
the sectional beam curvatures about the z- and y-ax.es, 
respectively [Bauchau and Craig (2009)]. The three un¬ 
knowns Si, K 2 , and K 3 are determined by the constraints 
that (i) the integral of over the midshaft cross 

section must give the total normal force N^, and (ii) the 
integral of the stress moment ( 0 , y, z) x must 

give the total bending moment M = M m (the axes 
origin in the {y, z) plane is set at the modulus-weighted 
centroid of the section, also called normal force cen¬ 
ter [Bauchau and Craig (2009)]). Explicit formulas for 
ei,K 2 , and K 3 as functions of N and M are pre¬ 

sented in Appendix C. We refer the reader to Bauchau 
and Craig (2009), Sec. 6.3 and Buenzli et al (2013) for 
their derivation. 

Determination of microscopic stress and strain distribu¬ 
tions based on micromechanical homogenisation theory 

Bone tissue stiffness C*'®®’^® is strongly influenced by 
the tissue’s microstructure, in particular its porosity, or 
equivalently, its bone volume fraction /bm- Bone volume 
fraction is a microstructural parameter defined at the 
tissue scale as the volume fraction of bone matrix in 
the RVE (Figure 1(d)): /bm = BV/TV = 1 — porosity, 
where BV is the volume of bone matrix in the RVE and 
TV is the tissue volume, i.e. the total volume of the RVE 
[Dempster et al (2013)]. In Buenzli et al (2013), we used 
an explicit power-law relationship C‘A®!r(/bm) (X fbJ 
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based on experimental relationships between bone stiff¬ 
ness and bone mineral content [Carter and Hayes (1977); 
Hernandez et al (2001)]. While regression approaches 
based on power-law relations are able to account for 
material properties in one principal direction, they are 
less accurate in estimating material properties in other 
principal directions. 

Here, we follow a different approach taken by Hellmich 
and colleagues using the framework of continuum mi¬ 
cromechanics [Hill (1963,1965); Zaoui (1997, 2002)]. Me¬ 
chanically, bone tissue can be considered as a two-phase 
material; a bone matrix phase (‘bm’) consisting of miner¬ 
alised bone matrix, and a vascular phase (‘vas’) consist¬ 
ing of vascular components, cells, extracellular matrix 
and other soft tissues present in Haversian canals and 
in the marrow. 

Continuum micromechanics provides a framework 
to estimate the tissue-scale stiffness tensor C‘‘®®“(/bm) 
from the microscopic stiffness properties of bone matrix 
and vascular pores, and assumptions on pore microar¬ 
chitecture and phase interactions [Hellmich et al (2008)]. 
The advantage of this approach is to provide (i) accurate 
three-dimensional estimates of and (ii) estimates 

of the microscopic stress and strain distributions of the 
bone matrix without recourse to costly micro-Hnite ele¬ 
ment analyses of the tissue microstructure [Fritsch et al 
(2009)]. Using the concept of continuum micromechan¬ 
ics is justihed in bone due to the separation of length 
scales between the RVE size and the characteristic sizes 
of the two-phase microstructures [Hellmich et al (2008); 
Scheiner et al (2013)]. We summarise below the premises 
upon which this approach is based. 

The tissue-scale stress and strain tensors and 

^tissue correspond to spatial averages over the RVE of 
the microscopic (cellular-scale) stress and strain tensors. 
Assuming that each phase within the RVE is homoge¬ 
neous, these spatial averages can be expressed as sums 
over the phases: 


_tissue/„ _ l \ _ ^ f _micro \ ' r _micro /oi 

o- <^y = 2^Jk(Tk , ( 3 ) 

j. V Jtv ^ 

et'®®“(r,t)= = (4) 


where fk{r,t) is the volume fraction of phase k (‘bm’, 
‘vas’), cr“‘®"'°(r,t) and e™‘®‘'°(r,t) are the microscopic 
stress and strain tensors in phase k. We emphasise that 
all these quantities still depend on the tissue-scale loca¬ 
tion r of the RVE in bone, whilst microscopic inhomo¬ 
geneities are encoded in the phase index k. It can be 
shown that due to the linearity of the constitutive equa¬ 
tions the phase strain tensor is related linearly 


with the tissue-scale strain tensor: 


^micro a . ^tissue 

= Ak : e 


( 5 ) 


where Ak is a fourth-order tensor called the strain con¬ 
centration tensor [Zaoui (2002); Hellmich et al (2008); 
Fritsch et al (2009)]. Assuming that Hooke’s law also 
holds for each phase at the microscopic scale, = 

c™'®™ : (with c™'®™ the stiffness tensor of phase 

k), one obtains from Eqs (3) and (5): 


_tissue 

CT 




^micro 


( \ '' £ ^micro .a \ , ^tissue _ 

2_^Jk‘^k : Afc J : e = 

k 


( 6 ) 


^tissue 


. ^tissue 


where 

C tissue £ ^micro .a , £ ^micro . a {'y\ 

— /bm ^bm ■ -^bm /vas ^vas * -^vas- (7) 

Equation (7) provides a relationship between the tissue- 
scale stiffness, C*'®®™, and the microscopic properties of 
the phases composing the tissue, /fc,®™'®™, and Ak- Be¬ 
cause mineral content across bone tissues only varies 
little on average [Scheiner et al (2013); Fritsch and 
Hellmich (2007)], c]]]])’"'® can be assumed constant and 
homogeneous, i.e., independent of r,t. The elastic mod¬ 
ulus (C™g®®° is likewise assumed independent of r,t and 
taken as that of water [Scheiner et al (2013)]. Both 
and c™*®®® have been measured experimentally, 
their values are listed in Table 1. The strain concentra¬ 
tion tensors Ak can be estimated by solving so-called 
matrix-inclusion problems of elasticity homogenisation 
theory, which use assumptions on phase shape within 
the RVE and phase interactions [Eshelby (1957); Laws 
(1977)]. For bone, accurate multi-scale homogenisation 
schemes were developed and validated experimentally 
[Hellmich et al (2008); Fritsch et al (2009); Morin and 
Hellmich (2014)]. These schemes provide explicit expres¬ 
sions for Ak depending on the phase volume fractions 
/bm and /vas- Because /vas = 1 - /bm, this defines both 
the /bm dependence of C**®®“® via Eq. (7), and a method 
to estimate the strains and stresses in the bone matrix 
at the microscopic level from those known at the tissue 
level: 

ebr°(r,t)=Abm(/bm) (8) 

^c®®(r,t) = cC“® : (Abm(/bm) : 

= Bbm(/bm) : <t‘‘®®®®, (9) 

where Hooke’s law (1) was used in the last equality 
in Eq. (9). The stiffness tensor C*‘®®®®(/bm), the strain 
concentration tensor Abm(/bm), and the stress concen¬ 
tration tensor ]Bbm(/bm) can be evaluated numerically 
at each location r in the femur midshaft cross section 
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and each time t based on the value of fhni{'r,t) and 
the expressions given in Fritsch and Hellmich (2007) 
and Scheiner et al (2013). 

Combined with beam theory, this procedure enables 
us to completely determine, at each time t, the spatial 
distribution across the femur midshaft of (i) the tissue- 
scale stress and strain tensors, and (ii) 

the microscopic stress and strain tensors of bone matrix, 

—micro —micro 
^bm ’ ^bm 

In this paper we will consider both the tissue-scale 
strain energy density (SED), and microscopic 

SED of the bone matrix, as local mechanical 

quantities sensed and transduced by bone cells. These 
SEDs are defined by: 

^tissue^^ _ J^^tissue . ^tissue . ^tissue ( 10 ) 

jTvmicro/^ ±\ _ 1 .micro . ^micro . .micro 

^bm • ^bm • ^bm • 

The SEDs defined in Eqs (10) and (11) will be used to 
formulate biomechanical regulation in the bone remod¬ 
elling equations. In the literature, biomechanical regula¬ 
tion is commonly based on the SED since this quantity 
is scalar and it integrates both microstructural state and 
loading environment [Fyhrie and Carter (1986); Mullen- 
der et al (1994); Ruimerman et al (2005)]. 


2.2 Bone tissue remodelling 

The tissue is assumed to be remodelled by a population 
of active osteoclasts (OCa) and active osteoblasts (OBa). 
Active osteoclasts are assumed to resorb bone at rate 
feres (volume of bone resorbed per cell per unit time). Ac¬ 
tive osteoblasts are assumed to secrete new bone matrix 
at rate fcform (volume of bone formed per cell per unit 
time). These cellular resorption and formation rates 
are taken to be constant and uniform. However, the 
bone volume fraction fhmifjt) of the tissue may evolve 
with site-specific rates depending on the balance be¬ 
tween the populations of active osteoclasts and active 
osteoblasts [Martin (1972); Buenzli et al (2013)]: 

^/bm(^5^) — feformOBa feresOCa- (12) 

In Eq. (12), OCa{r,t) and OBa(r,t) denote the average 
densities of active osteoclasts and active osteoblasts in 
the tissue located at r (number of cells in the RVE/TV, 
Figure 1(d)). The site-specific remodelling rate xbv(^, t) 
of the tissue at r (also called turnover rate) can be 
defined as the volume fraction of bone in the RVE that 
is resorbed and refilled in matched amount per unit time 
[Parfitt (1983), Sec. II.C.2.c.ii]. This corresponds to the 
minimum of the volume fraction of bone resorbed per 


unit time, fcresOCa, and volume fraction of bone formed 
per unit time, fcformOBa: 

XBv(^, t) = min{feresOCa, feformOBa}. (13) 

Any imbalance between resorption and formation in 
Eq. (12) is interpreted as surplus resorption or surplus 
formation with respect to the baseline of bone properly 
turned over in Eq. (13). 

Equation (12) enables us to track site-specific modi¬ 
fications of the midshaft tissue microstructure through 
fhmifjt), from which stress and strain distributions 
across the midshaft can be estimated at both the tissue 
scale and the microscopic, cellular scale, by means of 
Eqs (34)-(35) and (8)-(ll). 

2.3 Bone cell population model 

It remains to specify how the populations of active osteo¬ 
clasts OCa(r,t) and active osteoblasts OBa{r,t) evolve 
in the RVE located at r under mechanobiological, geo¬ 
metrical and biochemical regulations. Eor this, we use 
a continuum cell population model based on rate equa¬ 
tions, originally developed by Lemaire et al (2004), and 
later refined and extended by Pivonka and co-workers 
[Pivonka et al (2008, 2010); Buenzli et al (2012); Pivonka 
et al (2013); Scheiner et al (2013); Pivonka et al (2012)]. 

To highlight important biochemical couplings and 
regulations in osteoclastogenesis and osteoblastogene- 
sis, several differentiation stages of osteoclasts and os¬ 
teoblasts are considered. These biochemical interactions 
are mediated by several signalling molecules whose bind¬ 
ing kinetics are explicitly considered in the model, such 
as transforming growth factor |3 (TGFp), receptor-activator 
nuclear factor kB (RANK) and associated ligand RANKL, 
osteoprotegerin (OPG), and parathyroid hormone (PTH). 
The biochemical network of these couplings and regula¬ 
tions is summarised in Figure 3. 

Active osteoclasts (OCas) denote cells attached to 
the bone surface that actively resorb bone matrix. These 
cells are assumed to differentiate from a pool of osteo¬ 
clast precursor cells (OCpS) by the binding of RANKL to 
the RANK receptor, expressed on OCpS, which induces 
intracellular NFkB signalling. Osteoclast precursors are 
assumed to differentiate from a pool of uncommitted 
osteoclasts progenitors (OCu), such as haematopoietic 
stem cells, under the action of macrophage colony stim¬ 
ulating factor (MCSF) and RANKL signalling [Roodman 
(1999); Martin (2004)]. 

Active osteoblasts (OBas) denote cells at the bone 
surface that actively deposit new bone matrix. These 
cells are assumed to differentiate from a pool of os¬ 
teoblast precursor cells (OBps). This activation is in¬ 
hibited in the presence of TGFp. Osteoblast precursors 
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Fig. 3 Proposed cell population model of bone remodelling taking into account several developmental stages of osteoblasts 
and osteoclasts together with biochemical regulation, biomechanical regulation (via strain energy density, !?) and geometrical 
regulation (via the turnover function, XBv(/bm))- 


are assumed to differentiate from a pool of uncommit¬ 
ted osteoblasts progenitors (OBu), such as mesenchymal 
stem cells or bone marrow stromal cells, upon TGFp 
signalling [Roodman (1999)]. 

The rate equations governing the evolution of the 
tissue-average cell densities are given by: 

^OCp(r,t) =I?oc„ (mCSF, RANKL(tZ>', PTH)) OCu(/bm) 
- Pocp (RANKL(if', PTH)) OCp, (14) 

^OCa(r,t) =2?oCp (RANKL(tf^, PTH)) OCp 

-AloCp(TGFp)OCa, (15) 

^OBp(r,t) =I?oBp(TGFp)OBu(/b„0 +iPoBp(>Z') OBp 

-I?oBp(TGFp)OBp, (16) 

^OBa(r,t) =2?oBp(TGFp)OBp 

~ ^loBpOBa, (17) 

where Vi is the differentiation rate of cell type i (i = 

OCu, OCp, OBu, OBp) modulated by signalling molecules, 
Aoc^ is the apoptosis rate of active osteoclasts modu¬ 
lated by TGFp, is the (constant) apoptosis rate 

of active osteoblasts, Vob^ is the proliferation rate of 
osteoblast precursor cells, and ^ is the strain energy 
density, taken to be either or • 

The concentrations of the signalling molecules are 
governed by rate equations expressing mass action ki¬ 
netics of receptor-ligand binding reactions. Since time 
scales involved in cell differentiation and apoptosis are 
much longer than characteristic times of receptor-ligand 
binding reactions, the signalling molecule concentrations 
can be solved for in a quasi-steady state (adiabatic ap¬ 
proximation) [Buenzli et al (2012); Pivonkaet al (2012)]. 


Explicit expressions for the signalling molecules con¬ 
centrations and their modulation of the cell differentia¬ 
tion and apoptosis rates depending on receptor-ligand 
binding are presented in Appendix A.l and A.3. Below, 
we comment in more detail on new features of Eqs (14)- 
(17) that are included to model the geometrical and 
biomechanical feedbacks on bone cell populations. 

Geometrical feedback and turnover rate 

The local availability of bone surface to osteoclasts and 
osteoblasts is an important factor determining the propen¬ 
sity of initiating new bone remodelling events [Martin 
(1972); Buenzli et al (2013)]. A remarkable relationship 
between the density of bone surface, Sy, and bone vol¬ 
ume fraction, /bm, has been exhibited in bone tissues 
across wide ranges of porosities [Martin (1984); Fyhrie 
and Kimura (1999); Lerebours et al (2015)]. This prop¬ 
erty is particularly interesting from a computational 
modelling perspective as it enables to track microstruc- 
tural changes of bone tissues through the evolution of 
bone volume fraction only. 

In femur midshafts, bone tissue is usually compact, 
the bone volume fraction is high. However, during bone 
loss, bone volume fraction tends to decrease in the endo- 
cortical region. Due to the fact that /bm reaches values 
similar to trabecular bone volume fractions, this tissue 
has been called “trabecularised” cortical bone [Zebaze 
et al (2010)]. Here, we treat compact and porous tis¬ 
sues differently in terms of bone turnover rates [Parfitt 
(1983); Martin et al (1998)]. Different turnover rates in 
Eq. (13) can be achieved by assuming that OBu and 
OCu are functions of the bone matrix volume fraction. 
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Fig. 4 (a) Plot of the phenomenological relationship XBv(/bm) between turnover rate and bone volume fraction assumed in 
the model, the grey data points are the ones given by Parfitt [Parfitt (1983)]. (b) Dependence of OCu and OBu upon /bm 
assumed in the bone cell population model. 


/bin, which introduces a dependence of the active bone 
cell populations, OCa and OBa, upon /bm via Eqs. (14)- 
(17). This dependence may account both for the influ¬ 
ence of bone surface availability on turnover rate via the 
Sv(/bm) relation, and for influences of the biochemical 
microenvironment between cortical bone and trabecu- 
larised bone. 

Few experimental data explicitly associate turnover 
rate with microstructure. However, since the relation¬ 
ship Sv(/bm) is well established experimentally, it can be 
expected that a phenomenological relationship, XBv(/bm), 
associating bone turnover with bone volume fraction is 
well-defined. Parfitt reports that cortical bone of average 
bone volume fraction 0.95 has a turnover rate of 0.115 
cm^/day, corresponding to xbv(0-95) « 0.77 • 10“^/day 
with = 1.5 • 10® mm®. Moreover, he states that 

trabecular bone of average bone volume fraction 0.20 
has a turnover rate of 0.25 cm®/day, corresponding to 
Xbv(0.20) « 1.43 • 10-Vday with = 1.75-10® mm® 

[Parfitt (1983), Table 1 and Table 7]. We take for XBv(/bm) 
a dome-shaped function following Parfitt’s reported val¬ 
ues and having a zero turnover rate for /bm equal to 0 
and 1 (Figure 4(a)). The maximum of bone turnover is 
assumed to occur at /bm= 0.35, corresponding to typ¬ 
ical trabecular or trabecularised bone microstructures. 
These types of microstructures are expected to remodel 
at the highest rates due to the proximity of precursor 
cells in the marrow and the large availability of bone 
surface. 

The functions OCu(/bm) and OBu(/bm) are deter¬ 
mined such that the turnover rate obtained from the 
cell population model in a normal healthy state with 
balanced remodelling, matches the phenomenological re¬ 
lationship XBv(/bm). From Fqs (12) and (13), the bal¬ 


anced steady-state condition and the remodelling rate 
condition impose the constraints 

XBv(/bm) — ^formOPa(OCu(/bm), C)Pu(/bm)) (1^) 

= fcies^(OCu(/bm),OBu(/bm)) (19) 

at each value of /bm G [0,1], where the bar indicates 
steady-state values of the cell density variables. These 
two constraints were solved numerically with the turnover 
rate function XBv(/bm) reported in Figure 4(a) by using 
a trust-region dogleg method. The functions OCu(/bm) 
and OBu(/bm) obtained by this procedure are shown in 
Figure 4(b). These functions are used as input in Fqs 
(14)-(17) in all our simulations. This ensures that in 
steady state, each RVF of the midshaft cross section lo¬ 
cated at r remodels at rate XBv(/bm(»’)) without chang¬ 
ing its bone volume fraction. The functions OCu(/bm) 
and OBu(/bm) are assumed to hold unaffected in the 
various deregulating circumstances considered later on 
(i.e., osteoporosis and altered mechanical loading). 

The explicit calibration of the cell population model, 
Fqs (14)-(17), to site-specific tissue remodelling rates 
is a significant novelty compared to our previous tempo¬ 
ral model [Pivonka et al (2013); Scheiner et al (2013)]. 
This modification was made necessary to consistently 
describe the site-specific evolution of bone in the spatio- 
temporal framework of Buenzli et al (2013) whilst re¬ 
taining a cell population model that includes biochemi¬ 
cal regulations. 

Mechanical feedback and initial bone micro structure sta¬ 
bility 

A mechanical feedback is included in the cell population 
model such that underloaded regions of bone promote os- 
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teoclastogenesis and overloaded regions of bone promote 
osteoblastogenesis [Frost (1987, 2003)]. These responses 
are viewed as consequences of biochemical signals trans¬ 
duced from a mechanical stimulus sensed by osteocytes 
[Bonewald (2011)]. Osteocytes are known to express 
RANKL, which regulates osteoclast generation, and scle- 
rostin, which regulates osteoblast generation via Wnt 
signalling [Bonewald and Johnson (2008)]. Following 
Scheiner et al (2013), the resorptive response of the me¬ 
chanical feedback is assumed to act by an increase in the 
microenvironmental concentration of RANKL, whereas 
its formative response is assumed to act by an increase 
in the proliferation rate of osteoblast precursors. 

The exact nature of the mechanical stimulus sensed 
by osteocytes is still a matter of debate. It may include 
lacuno-canalicular extracellular fluid shear stress on the 
osteocyte cell membrane, extracellular fluid pressure, 
streaming potentials and direct deformations of the os¬ 
teocyte body induced by bone matrix strains [Knothe 
Tate (2003); Bonewald and Johnson (2008); Bonewald 
(2011)]. Due to the extensive network of osteocyte con¬ 
nections in bone [Buenzli and Sims (2015)], average bone 
matrix strains at a higher scale may also be sensed by 
the osteocyte network. Below, we assume that the me¬ 
chanical stimulus to the bone cell population model is 
described by a local strain energy density, <F(r,t). This 
strain energy density will be taken to be either the mi¬ 
croscopic, cellular-scale strain energy density of bone 
matrix, t), or the average tissue-scale strain en¬ 
ergy density, t), defined respectively in Eqs (11) 

and (10). 

To predict with our model the evolution of a real 
scan of midshaft femur under various deregulating cir¬ 
cumstances, it is important to assume that the bone 
scan represents a stable state initially in absence of any 
deregulation. In particular, this initial bone state is as¬ 
sumed mechanically optimal. This can be ensured by 
choosing the local mechanical stimulus acting onto the 
bone cells, /r(r,t), as a normalised difference between 
the current SED and the SED of the inital bone mi¬ 
crostructure 'F(r,0): 




'F(r, t) — !F(r, 0) 
W{r,0) + K 


( 20 ) 


The normalisation by 'F(r,0) in the denominator in 
Eq. (20) ensures that the mechanical stimulus is not 
over-emphasised away from the neutral axis where strain 
energy density takes high values. The small positive 
constant K = 1 ■ I0“® GPa is added to keep mechan¬ 
ical stimulus well defined near the neutral axis where 
>F(r, 0) « 0 (see also Discussion section 4). 


When negative, in Eq. (20) is assumed to 

promote the production rate of RANKL: 


/omech 

Prankl 




if/x(r,t)<0 
0 , if /x(r, t) > 0 


( 21 ) 


where «: is a parameter describing the strength of the 
biomechanical transduction (see section 2.6). This re¬ 
sults in increased RANKL signalling in underloaded con¬ 
ditions (see Eq. (28) in Appendix), and so in increased 
osteoclast generation in Eqs (14)-(15). 

When positive, /r(r, t) in Eq. (20) is assumed to pro¬ 
mote T^oBp, the proliferation rate of pre-osteoblasts in 
Eq. (16): 


{ 0 , if ^{r, t) <0 

^oBp • A • ^(r,t), if0<^(r,t)< 

^oBp, 

( 22 ) 

where A is a parameter describing the strength of the 
biomechanical transduction. The first term in (22) ac¬ 
counts for a transit-amplifying stage of osteoblast dif¬ 
ferentiation occurring in absence of mechanical stimu¬ 
lation [Buenzli et al (2012)]. The proliferation rate is 
assumed to saturate to the value "PcBp = 2PoBp in 
highly overloaded situations to ensure the stability of 
the population of OBpS [Buenzli et al (2012); Scheiner 
et al (2013)]. 

A similar type of mechanical feedback was imple¬ 
mented in purely temporal settings in Scheiner et al 
(2013). The initial strain energy density distribution, 
!F(r,0), is calculated from Eqs (lO)-(ll) and from the 
initial bone volume fraction distribution, /bm(»’,0), de¬ 
termined on the bone scan (described in the next sec¬ 
tion). 


2.4 Initial distribution of bone volume fraction from 
microradiographs 

The initial microstructural state of the midshaft bone 
cross section can be derived from high-resolution bone 
scans such as micro-computed tomography (microCT) 
scans or microradiographs. Since Haversian canals have 
an average diameter of about 40 /rm, at least 10 fxm reso¬ 
lution is required to evaluate intracortical bone volume 
fractions with sufficient accuracy. 

In the simulations presented in Section 3, we used the 
microradiograph represented in Figure 5(a) where the 
pixel size is 7 fxm. The femur sample was collected from 
a 21-year-old subject. The microradiograph was digi¬ 
tised and binarised by a thresholding operation based 
on pixel grey level. Bone matrix is assigned the value 1 
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(a) Microradiograph 
Posterior 



Anterior 


(b) Binarised image 




(c) Bone volume fraction 
distribution 

A 


-1 

I 0.9 

Jo.8 

0.7 
0.6 
0.5 
-0.4 
■ 0.3 



Fig. 5 (a) Microradiograph of a midshaft femur cross section (courtesy of C. David L. Thomas and John G. Clement, Melbourne 
Femur Collection), (b) Binarisation of the microradiograph and determination of the local /bm values: (i) at the periosteal 
region, (ii) in the intracortical region, and (iii) at the endosteal region; (c) Bone volume fraction distribution extracted from the 
radiograph and interpolated. The dashed line represents the location of the neutral axis. The origin of the coordinate system 
(j/, z) is taken at the normal force center, NC. The grey lines are the 10 mm along which we are studying the evolution of the 
model in the Results Section. 


irrespective of the degree of mineralisation, and intra¬ 
cortical pores are assigned the value 0. The distribution 
of the bone volume fraction, /bm(i’)O), across the mid¬ 
shaft was determined by calculating the volume of bone 
matrix in a disk of 2 mm diameter, centred at each pixel 
of the binarised image and divided by the disk’s area. 
For the points near the periosteal surface, only the por¬ 
tion of the disk contained into the subperiosteal area 
was used for this calculation (see Figure 5(b)). The dis¬ 
crete values of fhm defined at each pixel contained in 
the subperiosteal region were then interpolated into a 
continuous function, /bm(i',0), using Matlab’s 2D cubic 
interpolation procedure. The result is shown in Fig¬ 
ure 5(c). A similar exclusion was not performed at the 
endosteal surface since this surface is less well defined, 
in opposition to the periosteal surface, due to the pres¬ 
ence of ‘trabecular-like structures. Bone matrix volume 
fractions near the endosteal surface are averages of in¬ 
tracortical bone regions and regions in the bone marrow 
cavity. 


2.5 Numerical simulations 

The multiscale mechanobiological model of bone remod¬ 
elling presented in this paper is governed by a coupled 
system of (i) distributed ODEs describing the evolu¬ 
tion of bone cell populations at each location r in the 
midshaft femur (Eqs (14)-(17)); and (ii) non-local and 
tensorial algebraic equations determining the mechani¬ 
cal state of the tissue RVE at r, both at the tissue scale 


and at the microscopic scale (Eqs (2)-(ll)). The model 
is initialised with a bone volume fraction distribution 
across the midshaft femur deduced from high-resolution 
bone scans (Figure 5(a)) and with steady-state popu¬ 
lations of cells fulfilling the site-specific turnover rate 
conditions Eqs. (18)-(19). This initial state is thereby 
constructed to be a steady state of the model, in which 
the biochemical, geometrical and mechanobiological reg¬ 
ulations of resorption and formation are balanced. 


To solve this non-local spatio-temporal problem nu¬ 
merically, we use a staggered iteration scheme in which 
we first solve the mechanical problem (i.e., tissue-scale 
SED and microscopic SED) assuming constant mate¬ 
rial properties, and then solve the bone cell population 
model and evolve the bone volume fraction at each loca¬ 
tion r of the femur midshaft assuming constant mechani¬ 
cal feedback for a duration At. After Z\t, the mechanical 
problem is recalculated based on the updated bone vol¬ 
ume fraction distribution, fhmi'i', t-|-Z\t), and this proce¬ 
dure is iterated. The ODEs are solved using a standard 
stiff ODE solver (Matlab, ode 15s). The spatial discreti¬ 
sation is a regular grid with steps Ay = Az = 0.8 mm. 
Due to the separation of time scales between changes in 
the local mechanical environment (years) and changes 
in bone cell populations (days), the mechanical stimu¬ 
lus requires updating after durations Z\t = 2 years. The 
accuracy of the numerical results depending on Z\t is 
studied in Appendix B. 
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(a) Osteoporosis simulation via PTH infusion 



Time (year) 


(b) Spaceflight simulation 



Time (month) 


Fig. 6 (a) Evolution of the total bone mass in the cross section with time while simulating osteoporosis: calibration of the 
PTH infusion. Osteoporosis is characterised by a bone loss of 1%/year [Parfitt and Chir (1987); Nordin et al (1988); Szulc et al 
(2006)]. (b) Evolution of the total bone mass in the cross section with time while simulating a spaceflight mission: calibration 
of the mechanical parameters, A and k. 


2.6 Model calibration 

The model presented in this paper contains: (i) biome¬ 
chanical parameters associated with the estimation of 
<P'(r, t), and (ii) parameters associated with the bone cell 
population model. Biomechanical parameters as well as 
biochemical parameters were determined and validated 
in other studies [Scheiner et al (2013); Pivonka et al 
(2013, 2008); Buenzli et al (2012)] (Table 1). Here we 
calibrate the newly introduced parameters: (a) mechan¬ 
ical coupling parameters A and k (Eqs (21) and (22)), 
and (b) biochemical parameters related to the simula¬ 
tion of osteoporosis. 

Calibration of the hormonal deregulation for osteoporo¬ 
sis simulation 

In our previous temporal model [Scheiner et al (2013)], 
osteoporosis was modelled by an increase in systemic 
PTH together with a reduction in the biomechanical 
transduction parameters: A and k. In this paper, we 
simulate age-related bone loss using a single parameter 
perturbation, i.e., an increase in systemic PTH concen¬ 
tration. This increase is calibrated so as to obtain a loss 
of total bone cross-sectional area in the femur midshaft 
of 1% per year ^ [Parfitt and Chir (1987); Nordin et al 
(1988); Szulc et al (2006)]. The total bone cross-sectional 
area is defined by the integral of /bm(T, t) over the mid¬ 
shaft cross section. In the model, a rate of bone loss of 

^ The calibration is performed without mechanical adapta¬ 
tion (i.e. setting A = 0 and k = 0 in Eqs (21) and (22)) in order 
to compare both mechanical feedbacks in a more consistent 
way. 


1 %/year was obtained by an increase in systemic con¬ 
centration of PTH from 2.907 pM to 2.954 pM (1.62% 
increase) (see Figure 6(a)). 

Calibration of mechanobiological feedback 

The rate of change in bone mass due to mechanical 
feedback is determined in the model by the biomechan¬ 
ical transduction parameters A (in Eq. (22)) and n (in 
Eq. (21)). To calibrate these parameters, we used data 
gathered from mechanical disuse and re-use experiments. 
It has been shown that cosmonauts undertaking long 
mission space flights lose bone mass at a rate of ap¬ 
proximately 0.3% per month [Vico et al (2000)]. This 
microgravity-induced bone loss is only slowly recovered 
after return to Earth. No significant bone gain is ob¬ 
served after 6 month exposure to normal gravity on 
Earth [Vico et al (2000); Collet et al (1997)]. 

In our multiscale model, microgravity is simulated 
as a 80% reduction of the normal mechanical loads ex¬ 
perienced by the femur, i.e., TVmicrogravity = 0.2Af and 
ATmicrogravity = 0.27VT. Based on these reduced loads, 
the parameter n has been calibrated such that 1.8% of 
total bone cross-sectional area is lost after 6 months. 
We found such rate of loss with «: = 18 pM/day when 
the mechanical stimulus is based on the microscopic 
SED, tf'^“°(r,t) (see Figure 6(b)), and k = 19 pM/day 
when the mechanical stimulus is based on the tissue- 
scale SED, t). After return to Earth, rates of 

bone recovery are too low to be detected after 6 months 
[Collet et al (1997)]. We performed a parametric study 
investigating various strengths of A. Using parameter 
values of A > 1 in our model results in significant bone 
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gain after 6 months, while A < 1 results in small bone 
gain. Based on these results we use A = 0.5 for both the 
microscopic and tissue-scale mechanical stimuli. 

3 Results 

In this section we present numerical simulations of the 
evolution of the midshaft femur cross section subjected 
to either: (i) changes in mechanical environment (Sec¬ 
tion 3.1), or (ii) hormonal deregulation simulating os¬ 
teoporosis (Section 3.2). We also investigate how site- 
specific bone loss may depend on whether mechanical 
stimulus is sensed at the microscopic, cellular scale, or 
at the tissue scale. 

3.1 Bone loss due to mechanical disuse 

Figure 7(a) represents site-specific changes of the femur 
midshaft cross section simulated by the model assuming 
a 80% reduction in the normal mechanical loading. This 
reduction in mechanical loading may represent micro¬ 
gravity in long spaceflight missions (see Section 2.6) or 
prolonged bed rest. The Figure depicts the difference 
between the bone volume fraction distribution after 6 
months of mechanical disuse and the initial bone vol¬ 
ume fraction. It can be seen that bone loss is site-specific 
with more bone loss occurring near the endosteal sur¬ 
face. Close to the neutral axis, only limited loss of bone 
is observed. 

3.2 Simulation of osteoporosis due to hormonal 
deregulation 

Figures 7(b) and (c) represent the site-specific changes 
of the midshaft cross section that occur after 40 years 
of simulated osteoporosis when the mechanical feedback 
acting onto the bone cell population model is based on 
the microscopic SED. Figure 7(b) depicts the difference 
between the /bm distribution at the end of the simu¬ 
lation and the initial distribution. Figure 7(c) depicts 
the /bm distribution at the end of the simulation. Bone 
loss occurs everywhere in the cross section except at the 
medial and lateral sides. The loss is site-specific with 
higher rates of loss in the endocortical region and around 
the neutral axis, close to the antero-posterior axis. This 
pattern of bone loss is consistent with the high poros¬ 
ity commonly observed in these regions in osteoporotic 
subjects (Figure 7(d), arrows). The simulation exhibits 
a sharp transition between a very porous endocortical 
region and a dense intracortical region towards the pe¬ 
riosteum. Although perhaps less pronounced, such a 


transition is also observed in the microradiograph of 
Figure 7(d) (dashed lines). In contrast to the osteo¬ 
porosis simulation, the simulation of mechanical disuse 
(Figure 7(a)) shows that bone was lost all over the cross 
section, with little change around the neutral axis. In 
both simulations, bone was lost predominantly in the 
endocortical region. 

In Figure 8, we show how the distributions of the 
following quantities evolved along the y- and the z-axes 
during the simulation of osteoporosis: (a) the bone vol¬ 
ume fraction, (b) the microscopic mechanical stimulus, 
used as mechanical stimulus in this simulation, 
(c) the tissue-scale mechanical stimulus, not used 

as mechanical stimulus in this simulation, and (d) the 
turnover rate. Along both axes, the regions in which 
bone volume fraction transitions from low to high (3 
< y < 6 mm and -7 < z < -5 mm) are resorbed at 
higher rate, due to the higher values of xbv in these 
regions (Figure 8(d)). As time progresses, bone volume 
fraction is strongly reduced in the endocortical region, 
leading to an expansion of the marrow space and a re¬ 
duction in cortical wall width. This is accompanied by 
a shift of the maximum of Xbv towards the periosteum. 
Along the y-axis (near the neutral axis), bone is lost 
at a high rate not only in the endocortical region but 
also near the periosteum, as can be seen by the gradual 
increase in turnover rate in the whole cortical width 
(Figure 8(d)). In contrast, along the z-axis, bone is lost 
at a high rate only at the endosteum where turnover 
rate maintains a well-defined peak. The intracortical 
region (z < -7 mm) is preserved even after 40 years of 
simulated osteoporosis. 

Microscopic vs tissue-scale mechanical stimulus 

Comparing Figures 8(b) and (c), we can observe that the 
values of the mechanical adaptation stimuli are strongly 
dependent on the length scale at which they are calcu¬ 
lated, i.e. tissue scale or microscopic scale. Along the 
z-axis, is always positive (Figure 8(b)), whereas, 

^tissue negative values in the endocortical region 

(Figure 8(c)). Regions with high values of and 

positive values of correlate with regions where 

the bone matrix is preserved. Regions with low values 
of and negative values of correlate with 

regions where the bone matrix is resorbed. The Figures 
also show a qualitative and quantitative difference in 
mechanical stimuli and between the y- and 

z-axes. The mechanical stimulus is asymmetric between 
the antero-posterior axis and lateral-medial axis due to 
the assumed bending loading state. Along the y-axis, 
no important variation can be observed between endo¬ 
cortical and periosteal regions. Along the z-axis, both 
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(a) After spaceflight:- 6m) -/bmCl ^ 0) (b) After simulated osteoporosis: 

/bm(f = 40y) -./bm(f = 0) 




(d) Microradiograph of a 89 yo individual 


-5 0 5 

z - axis (mm) 


(c) After 40 years of simulated osteoporosis 


0.6 

0.4 



Fig. 7 (a) Difference between the bone volume fraction distribution in the cross section after a 6-months spaceflight mission 
with the initial distribution, (b) Difference between the bone volume fraction distribution in the cross section after a 40-years of 
simulated osteoporosis with the initial distribution, (c) Bone Volume Fraction distribution in the cross section after 40 years of 
simulated osteoporosis, (d) Microradiograph of a human femur cross section from an 89 years old individual. The dashed lines 
highlight regions with sharp transition between porous and compact tissue. The arrows point out regions with high porosity 
along the antero-posterior axis. 



stimuli exhibit much lower values in the endocortical 
region than at the periosteum or in the marrow cavity. 
We note that the mechanical stimulu are not zero in 
the marrow even when /bm = 0 due to the assumed 
vascular stiffness. 

Figure 9 compares the evolution of bone volume 
fraction during simulated osteoporosis when the me¬ 
chanical stimulus acting onto the bone cells is either (i) 
absent, (ii) based on the microscopic mechanical stimu¬ 
lus, or (iii) based on the tissue-scale mechanical 

stimulus, All cases exhibit strong endocortical 

bone loss with little difference in the expansion rate of 
the medullary cavity. A slightly steeper endosteal wall 
is created along the z-axis during the simulation using 

3 For the simulation in case (iii), the mechanical transduc¬ 
tion strength parameters are: k = 19 pM/day and A = 0.5, 
calibrated with the tissue-scale mechanical stimulus while sim¬ 
ulating spaceflight. 


tissue-scale mechanical stimulus, and a region with very 
low bone volume fraction (/bm — 0.1) is preserved in 
the medullary cavity during the simulation with micro¬ 
scopic mechanical stimulus. Intracortical bone towards 
the periosteum is preserved along the z-axis by both me¬ 
chanical stimuli, but it is resorbed more strongly along 
the y-axis in the simulation with tissue-scale mechanical 
stimulus. 


4 Discussion 

Endocortical hone loss 

The loss of endocortical bone, with its associated ex¬ 
pansion of the marrow cavity and cortical wall thinning, 
is a trait common to several bone disorders and dereg¬ 
ulations of remodelling. It is observed in osteoporosis 
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Fig. 8 Evolution of (a) Bone volume fraction; (b) Mechanical stimulus, /j.(r, t), at the microscopic scale; (c) Mechanical stimulus, 
at the tissue scale; and (d) Turnover rate along the y and z-axis during the simulation of osteoporosis. 


[Feik et al (1997); Parfitt (1998); Bousson et al (2001); 
Thomas et al (2005); Szulc et al (2006); Zebaze et al 
(2010)], vitamin D deficiency [Busse et al (2013)], hyper¬ 
parathyroidism [Hirano et al (2000); Burr et al (2001); 


Turner et al (2002)], but also during disruptions of nor¬ 
mal mechanical loading, such as in prolonged bed rest 
[Leblanc et al (2007); Rittweger et al (2009)], long term 
space missions [Vico et al (2000); Lang et al (2004)], 
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— Initial distribution 

— Without mechanical feedback 

— Microscopic scale mechanical stimulus 

— - - Tissue-scale mechanical stimulus 



Fig. 9 Bone volume fraction along the y- and 2 :-axes, the initial distribution and the distributions after 40 years of simulated 
osteoporosis, without mechanical regulation, with mechanical regulation based on microscopic SED, and with mechanical 
regulation based on tissue-scale SED. 


trauma-induced paralysis such as spinal cord injury [Ki- 
ratli et al (2000); Eser et al (2004)], and as well in ani¬ 
mal studies: muscle paralysis [Warner et al (2006); Ausk 
et al (2012, 2013)] or hind-limb disuse induced by tail 
suspension [Bloomfield et al (2002); Judex et al (2004)]. 

Our numerical simulations of osteoporosis and me¬ 
chanical disuse are consistent with these experimental 
findings. All Figures 7, 8 and 9 highlight the strong 
site-specificity of bone loss under deregulations of bone 
remodelling. The endocortical region systematically un¬ 
dergoes the most significant loss. This similarity arises 
despite the fact that in the simulation of mechanical 
disuse, the deregulation is non-uniform in the cross sec¬ 
tion (due to the uneven distribution of mechanical loads) 
whereas in the simulation of osteoporosis, the hormonal 
deregulation is uniform in the cross section. 

The precise mechanisms that underlie the predom¬ 
inant loss of bone in the endocortical region are still 
poorly understood [Raisz and Seeman (2001); Thomas 
et al (2005); Squire et al (2008); Ausk et al (2012)]. Me¬ 
chanical adaptation has been suggested as a potential 
mechanism [Frost (1997); Burr (1997); Thomas et al 
(2005); Jepsen and Andarawis-Puri (2012)]. Bone loss 
induced by mechanical disuse redistributes mechanical 
loads towards the periosteum, where bone volume frac¬ 
tion is higher. This could unload endocortical regions 
and thereby accelerate their resorption. Reduced physi¬ 
cal activity and muscle strength in ageing subjects sup¬ 
port this hypothesis [Frost (1997)]. However, the ubiq¬ 
uity of endocortical bone loss in situations in which me¬ 
chanical loading is not significantly modified suggests 
that other mechanisms are at play. The morphological 
influence of the tissue microstructure on the rate of 
bone loss has been suggested to be another important 


factor [Martin (1972); Squire et al (2008); Zebaze et al 
(2010); Buenzli et al (2013)]. Cortical bone has little 
bone surface available to bone cells, but this surface 
expands during bone loss, which could increase the ac¬ 
tivation frequency of remodelling events. If remodelling 
is imbalanced, this may lead to an acceleration of bone 
loss, and to an increase of the available surface until the 
tissue microstructure becomes so porous that its surface 
area reduces with further loss [Martin (1972); Raisz and 
Seeman (2001)]. 

We have shown previously the possibility of this 
morphological mechanism to explain cortical bone tra- 
becularisation in both temporal [Pivonka et al (2013)] 
and spatio-temporal settings [Buenzli et al (2013)]. The 
spatio-temporal model proposed in the present work 
incorporates both mechanical adaptation and a mor¬ 
phological feedback of the microstructure on turnover 
rate. In Figure 9, our simulations of osteoporosis con¬ 
ducted with and without mechanical feedback suggest 
that the rate at which the medullary cavity expands 
and the cortical wall thins is only marginally dependent 
on mechanical adaptation. This rate is primarily due 
to the high turnover rates present in the endocortical 
region (Figure 8(d)), i.e., due to the morphological influ¬ 
ence of microstructure on the rate of loss. This proposed 
mechanism is consistent with the observation that dis¬ 
tinct conditions exhibit endocortical bone loss, whether 
mechanical loading is disrupted or not. 

Model formulation of morphological feedback 

In the cell population model of Pivonka et al (2013), 
the morphological influence of the tissue microstructure 
was included through the specific surface of the tissue 
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[Martin (1984); Lerebours et al (2015)] normalised by 
its initial value. This normalisation allowed to maintain 
the same cell behaviour in both cortical and trabecular 
bones. However, it leads to a turnover rate that is ini¬ 
tially independent of bone volume fraction, and so the 
same in cortical and trabecular bone. The morpholog¬ 
ical feedback proposed in the present model differs by 
(i) avoiding a dependence on an initial reference state 
(i.e., absence of normalisation to allow microstructure- 
dependent turnover rates), and (ii) by referring to turnover 
rate (a dynamic biological quantity) instead of specific 
surface (a morphological characterisation of the microstruc¬ 
ture). 

Whilst specific surface can be estimated directly 
from high-resolution scans of bone tissues [Chappard 
et al (2005); Squire et al (2008); Lerebours et al (2015)], 
quantitative links between Sy and cell numbers remain 
unclear [Martin (1972); Parfitt (1983); Pivonka et al 
(2013)]. The direct reference to turnover rate, in the 
present model, makes the model more accurate, due 
to the straightforward link between turnover rate and 
cell populations (see Eq. (13)). Unfortunately, turnover 
rate is rarely measured experimentally by cell counts 
or volumes of bone resorbed and re-formed. It is most 
commonly characterised by measurements of serum con¬ 
centrations of bone resorption and/or formation mark¬ 
ers [Szulc et al (2006); Burghardt et al (2010); Mal- 
luche et al (2012)], which are difficult to relate quan¬ 
titatively to cell numbers or bone volume at a partic¬ 
ular bone site. Whilst the phenomenological relation¬ 
ship XBv(/bm) that we assumed between turnover rate 
and microstructure remains to be studied quantitatively, 
such a relationship has been suggested in several studies 
[Felsenberg and Boonen (2005); Burghardt et al (2010); 
Malluche et al (2012)]. 

Nature of the mechanical stimulus 

The nature of the mechanical stimulus sensed by bone 
cells and transduced into signals prompting resorption 
or formation has been a matter of discussion for many 
years. A number of computational studies simulating me¬ 
chanical adaptation of bone microstructure suggested 
that the strain energy density could be a good candi¬ 
date. Ruimerman et al (2005) tested several mechanical 
stimuli and concluded the SED gave best results when 
comparing simulations outcomes with biological param¬ 
eters such as porosity, trabecular number or adaptabil¬ 
ity to external loading. However, Levenston and Carter 
(1998) argued that the drawback of using the SED is 
that it does not lead to a different response when bone 
is stimulated in tension or in compression. In the litera¬ 
ture, most computational models use the SED because 


it is a scalar representing both microstructure and me¬ 
chanical loading [Fyhrie and Carter (1986); Mullender 
et al (1994); Van Rietbergen et al (1999); Van Oers 
et al (2008); Scheiner et al (2013)]. Quantitative criteria 
based on experimental observations are still lacking, es¬ 
pecially ones testing the tensorial aspects of mechanical 
loading conditions. For our purpose of studying tissue- 
scale average changes in bone volume fraction, these 
tensorial aspects are likely to be secondary. Hence, we 
have based our mechanical stimulus on the strain en¬ 
ergy density (see below for a discussion of the scale). 
We note here that other mechanical quantities have 
also been proposed and studied for their magnitude and 
possible influence onto osteocytes, such as fluid shear 
stress and fluid pressure in the lacuno-canalicular sys¬ 
tem [Knothe Tate et al (1998); Burger and Klein-Nulend 
(2003); Tan et al (2007); Bonewald and Johnson (2008); 
Adachi et al (2009b)]. 

Mechanical adaptation also relies on the comparison 
of the current mechanical state with a reference state. 
The definition of this reference state remains unclear 
[Frost (1987); Carter and Beaupre (2001)]. Our choice 
is to take as mechanical reference state the initial dis¬ 
tribution of the strain energy density in the midshaft 
femur. This choice introduces a memory of the stimu¬ 
lus “normally” experienced in a certain region of the 
tissue. This memory effect leads to a position-dependent 
reference state which can be interpreted as taking into 
account different sensitivities of the mechano-sensing 
cells depending on where they are located [Skerry et al 
(1988); Turner et al (2002); Robling et al (2006)]. 

Neutral axis and site-specific hone adaptation 

A common issue in models of mechanical adaptation is 
the risk to resorb too much bone in regions that are nat¬ 
urally unloaded. Such regions may exist when bending 
moment is large enough with respect to compressive or 
tensile forces. In the human midshaft femur, a neutral 
axis runs approximately along the antero-posterior axis 
[Lanyon and Rubin (1984); Gordey and Gautier (1999); 
Thomas et al (2005); Martelli et al (2014)]. To prevent 
excessive resorption in such regions, some models have 
considered torsional loads [Van der Meulen et al (1993); 
Garpenter and Garter (2008)], average values of periodic 
dynamic loads under which the neutral axis moves ]Van 
der Meulen et al (1993); Garter and Beaupre (2001)], 
or a residual background of mechanical stimulus mod¬ 
elling muscle twitching and other background mechani¬ 
cal forces [Mittlmeier et al (1994); Garpenter and Garter 
(2008)]. 

Such additional features were not introduced explic¬ 
itly in our model. The strength of the mechanical stim- 
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ulus around the neutral axis remained weak in our sim¬ 
ulation of mechanical disuse (Figure 7(a)). This is due 
to the fact that stimulus sensitivity is prescribed ac¬ 
cording to the initial state. The neutral axis did not 
move substantially during the simulation, and so the 
difference in strain energy density remained small. Re¬ 
sorption around the neutral axis was pronounced in our 
simulation of osteoporosis (Figures 7(b), (c) and 8) due 
to hormone-induced remodelling imbalance. Resorption 
was limited by the duration of the simulated osteoporo¬ 
sis (40 years) and the calibration of the overall bone loss 
according to experimental data. 

The bending moment exerted onto the femur at the 
midshaft creates a strong asymmetry in the local me¬ 
chanical state. Over time, this asymmetry leads to a 
cortical wall thickness which differs between the j/-axis 
(antero-posterior axis) and the z-axis (lateral-medial 
axis), as seen in Figure 8(a). Asymmetries in cortical 
wall thickness and bone volume fraction, in osteoporotic 
patients, are commonly observed (Figure 7(d)) [Feik 
et al (2000); Thomas et al (2005); Zebaze et al (2010)]. 

Microscopic vs tissue-scale mechanical regulation 

Mechanical deformations of bone matrix can be sensed 
by osteocytes at the microscopic, cellular scale by defor¬ 
mation of the cell body, transmitted either through di¬ 
rect contact with the matrix, or through changes in fluid 
flow or hydrostatic pressure [Weinbaum et al (1994); 
Knothe Tate (2003); Adachi et al (2009b,a); Bonewald 
(2011)]. However, osteocytes are highly connected to 
one another and to other cells in the vascular phase 
by an extensive network [Marotti (2000); Kerschnitzki 
et al (2013); Buenzli and Sims (2015)]. Whilst signal 
transmission mechanisms in this network remain to be 
determined, it is possible that the network integrates de¬ 
formations of both the matrix and vascular phases before 
transducing them into a biochemical response, enabling 
a mechanical sensitivity of the network to tissue-average 
stresses and strains. 

The uncertainty of the scale at which mechanical 
stimulus is sensed in bone has motivated many compu¬ 
tational studies to estimate stress concentration effects 
in bone microstructures [Hipp et al (1990); Kasiri and 
Taylor (2008); Gitman et al (2010)]. However, few stud¬ 
ies have explored the changes that occur during simu¬ 
lated bone loss when using microscopic or tissue-scale 
mechanical stimulus. 

Our simulations of osteoporosis show that most of 
the difference between the mechanical stimulus at the 
microscopic and tissue scales occurs near the endosteum 
and neutral axis (Figure 8(b,c)). Changes in bone vol¬ 
ume fraction were similar in both simulations. Stress 


concentration effects captured in the microscopic me¬ 
chanical stimulus (but not in the tissue-scale mechan¬ 
ical stimulus) resulted in maintaining a region of low 
porosity (/bm — 0.1) near the medullary cavity and in 
widening the transition between endocortical and intra- 
cortical bone volume fractions (Figure 9). 

Osteoporotic human femur midshafts exhibit a wide 
range of variability, reflecting the multiple factors in¬ 
fluencing bone loss [Feik et al (1997, 2000); Thomas 
et al (2005); Zebaze et al (2010)]. The expansion of 
the medullary cavity and thinning of the cortical wall 
are commonly reported, but other changes in midshaft 
tissue microstructures have been studied less system¬ 
atically. Depending on the subject and their specific 
condition, the transition between porous endocortical 
bone and dense intracortical bone may be sharp or wide, 
and highly porous microstructures near the endosteum 
may be found or not [Feik et al (1997), Figure 6]. 

Our model possesses several limitations which pre¬ 
vent at this stage to draw definite conclusions about 
the mechanical regulation of the tissue. The mechanical 
state is calculated only based on bone volume fraction. 
Other microstructural parameters such as the connec¬ 
tivity of the microstructure are not accounted for. Loss 
of connectivity is observed in osteoporotic trabecular 
bone [Parfitt et al (1987); Mosekilde (1990); Raisz and 
Seeman (2001)], which could lead to mechanical disuse 
and so increase in resorption. Periosteal apposition is of¬ 
ten reported and believed to result from a compensation 
of endocortical bone loss in osteoporotic patients [Szulc 
et al (2006); Russo et al (2006); Jepsen and Andarawis- 
Puri (2012)]. Our simulations assumed the periosteal 
surface to be fixed, which could limit the expansion rate 
of the medullary cavity. Finally, our simulation of os¬ 
teoporosis assumed a constant level of physical activity. 
A reduction in physical activity with age could further 
limit the preservation of bone matrix by mechanical 
feedback. 


Summary and conclusions 

In this paper a novel spatio-temporal multiscale model 
of bone remodelling is proposed. This model bridges 
organ, tissue and cellular scales. It takes into account 
biochemical, geometrical, and biomechanical feedbacks. 
The model is applied to simulate the evolution of a 
human femur midshaft scan under mechanical disuse 
and osteoporosis. It enables us to investigate how these 
scales and feedbacks interact during bone loss. Our nu¬ 
merical simulations revealed the following findings: 
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— Endocortical bone loss during both osteoporosis and 
mechanical disuse is driven to a large extent by site- 
specific turnover rates. 

— Mechanical regulation does not influence significantly 
the expansion rate of the medullary cavity. 

— Mechanical regulation helps preserve cortical bone 
near the periosteum. It explains site-specific differ¬ 
ences in the bone volume fraction distribution in the 
midshaft cross section during osteoporosis such as 
increased porosity near the neutral axis, and thicker 
cortical wall along the medial-lateral axis of the fe¬ 
mur midshaft, due to the anisotropy of the mechan¬ 
ical stimulus in the presence of bending moments. 

— The inclusion of stress concentration effects in the 
mechanical stimulus sensed by the bone cells has 
a pronounced effect on porosity in the endocortical 
region. 


Our methodology provides a framework for the future 
development of patient-specific models to predict loss 
of bone with age or deregulating circumstances. 


A Complements on the model description 


A.l Differentiation rates and signalling molecules in 
the cell populations model 

In Section 2.3, we presented the simplified equations of the 
bone cells population model. Here are the developments of 
these equations. 
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In those equations, several signalling molecules play a role: 
TGFp, RANK, RANKL, OPG, MCSF and PTH. The concentra¬ 
tions of these molecules follow the principles of mass action 
kinetics of receptor-ligand reactions. Due to the separation 
of scale between the cells differentiation and apoptosis rates 
and the receptor-ligand binding reactions, we solve them in a 
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A.2 Parameter values 

See Table 1. 


A.3 Recalibration of the model 


Since OGu and OBu vary with /bm so as to retrieve experi¬ 
mentally valid turnover rates, some other parameters required 
modification compared with previous versions of the cell pop¬ 
ulation model in which OCu and OBu were constant and un¬ 
calibrated [Buenzli et al (2012); Pivonka et al (2013)]. 

By comparing the cell densities between this model and 
the previously published one [Buenzli et al (2012)], we can 
determine scaling coefficients which allows a systematic cali¬ 
bration of ) and ( gitSitg ) ■ Indeed, these func- 

^OCa 

tions depend on the active and precursor cell densities. In the 
original models, the constants in these functions were cali¬ 
brated such as to obtain a strong biochemical feedback re¬ 
sponse. Maintaining this strong biochemical response is the 
aim of this re-calibration. 

The calibration is realised at /bm= 0.90. Both the turnover 
rate value and the values of fcres and fcform) have been changed 
according to the literature. Hence, by isolating OBa and OCa 
in the two new constraints of the steady state, the active 
osteoblast and active osteoclast read: 
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if (5 is the coefficient of proportionality between the new bone 
turnover rate and the previous one; /3 = S ■ kres/k^^^ and 
7 = (S-fcform/fcfoyjii' These coefficients are introduced in the de¬ 
termination of TGF/3 and OPG. Previously, TGF/3 was [Buen- 
zli et al (2012)]: 
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The new one becomes: 
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Table 1 Nomenclature 


Symbol Description 

Value 


Xbv Turnover rate 

OBu Uncommitted osteoblasts 

OCu Uncommitted osteoclasts 

^form Daily volume of bone matrix formed per osteoblast 

kres Daily volume of bone matrix resorbed per osteoclast 

A Strength of the mechanical transduction in formation 

K Strength of the mechanical transduction in resorption 

Stiffness tensor of the bone matrix phase 


Stiffness tensor of the vascular phase 


M 


Normal force 
Bending moment 


Extrapolated function of Am from [Parfitt (1983)] 

Given function of fbm, determined to fulfil the steady state 
Given function of Ami determined to fulfil the steady state 
150 /rm^/day [Buenzli et al (2014)] 

9.43T0^/im^/day [Buenzli et al (2014)] 

0.5 (Parametric study) 

18 pM/day (with 19 pM/day (with (Parametric study) 

/ 28.4 11.0 10.4 0 0 0 


0 

11.0 20.8 10.3 0 
10.4 10.3 18.5 0 
0 0 0 12.9 

0 0 0 0 
0 0 0 0 
/I 1 1 0 0 ON^ 


11.5 0 
0 9.3 J 


GPa [Ashman et al (1984); Fritsch and Hellmich (2007)]^ 


2.3 ■ 


1110 0 0 
1110 0 0 
000000 
000000 
Vo 0 0 0 0 o/ 


GPa [Murdock (1996)] 


-700 N (see Section 2.1) 
50 Nm (see Section 2.1) 


DoCu 

Differentiation rate of OCu 

into OC] 

Dob,, 

Differentiation rate of OBu 

into OB] 

DoCp 

Differentiation rate of OCp 

into OC; 

DoBp 

Differentiation rate of OBp 

into OB; 

PoBp 

Proliferation term of OBp 


^OCa 

Apoptosis rate of OCa 


^OBa 

Apoptosis rate of OBa 



0.42/day [Pivonka et al (2013)] 
0.7/day ^ 

2.1/day 

0.166/day 

3.5-10-3/day 

5.65/day 

0.211/day 


^TGF^ Parameter for TGF(3 binding on OBu and OCa 
Parameter for TGFp binding on OBp 
Parameter for PTH binding on OB (activator) 
^OB^ Parameter for PTH binding on OB (repressor) 

Parameter for RANKL binding on OCp 
Number of RANK receptors per OCp 
^oc Parameter for MCSF binding on OCu 

PpTH Systemic concentration of PTH 

PpTH PpTH when simulated osteoporosis 

Dqpg Degradation rate of OPG 

/^OB*^ Production rate of OPG per OBa 

OPGsat Saturation of OPG 

Drankl Degradation rate of RANKL 

Production rate of RANKL per OBp 
^OPG^^ Parameter for RANKL binding on OPG 
^RANK^ Parameter for RANKL binding on RANK 
^TGF^ Degradation rate of TGFp 

Density of TGFp stored in the bone matrix 


5.6310-^ pM 
1.89-10-3 pM 
150 pM 
0.222 pM 
16.65 pM 
1 - 10 ^ 

1- 10-3 pM [Pivonka et al (2013)] 
2.907 pM 

2.954 pM 
0.35/day 
1.63-10®/day 

2- 103 pM 
10/day 
1.68-103/day 
1-10-3/pM 
0.034/pM 

2/day 
1-10-2 pM 


^ Note that in comparison with Scheiner et al (2013), the x- and 2 :-axes are switched. 
3 Unless otherwise specified, parameter values are taken from [Buenzli et al (2012)] 
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The same manipulation is realised on the determination of C Generalised Beam theory for inhomogeneous 
OPG. The factor/3gPG OBa in Eq. (27), becomes OB!;®"'7-ijnaterials 


B Update frequency of mechanical state in the 
numerical algorithm 

In our model, to simulate osteoporosis and the change of poros¬ 
ity with time, we need to solve the temporal equations of 
the bone cell populations model, Eqs (14)-(17) and Eq. (12). 
Those equations via the mechanical feedback are correlated 
to the spatial Equations (2) and (42). Knowing the porosity 
distribution is required to determine the stress and strain dis¬ 
tributions. Hence we have a semi-coupled algorithm (Figure 
2 ). 

However, due to the separation of time scale we can de¬ 
compose the problem into two parts. Indeed, it takes more 
time for the microstructure to change significantly enough to 
influence the bone cell populations model. Therefore, we solve 
the bone cell populations model for a duration At, assuming 
the mechanical feedback to be constant in this time interval. 
Then, we recalculate the stress and strain distributions based 
on the new porosity distribution, and this becomes the new 
mechanical feedback. 

A sensitivity analysis of the solution in the time step At 
of evolution of cell densities and bone matrix volume fraction 
is required. For very small time steps (At < 1 day) one would 
expect the algorithm to converge to the exact solution. On 
the other hand for very large time steps (At > 5 years) a 
large deviation from the exact solution is expected. Figure 10 
shows the evolution of the bone matrix volume fraction for 
one selected RVE [y = 0, z = -10 mm) in the cross section. 
These simulations show that time steps of At = 250 days, 1 
year and 2 years lead to very similar evolution of the bone 
matrix volume fraction. On the other hand. At = 5 years 
and 10 years lead to strong deviations from the smaller time 
increments. For all the simulations of 40 years of osteoporosis, 
we used a time step of 2 years. 



Fig. 10 The evolution of the bone matrix volume fraction for 
different time steps. Note that this RVE is in the intracortical 
region which undergoes first resorption then formation due to 
the redistribution of the mechanical loads. 


In the following, we represent the governing equations using 
a Cartesian {x, y, z) coordinate system. The a:-axis represents 
the beam axis and coincides with the direction of the vascular 
pores (i.e., Haversian systems). The y and 2 coordinates de¬ 
scribe a material point in the cross section (Figure 1(c)). The 
origin of the system is known as the Normal force Center: NC. 
Since our cross sections are inhomogeneous all the quantities, 
including the stiffness, are functions of y and 2 . 

First, based on the constitutive relation: Hooke’s law, we 
determine the strain and stress relation: 

<T*'==“®(j/, 2 ,t) = C*'==“®(j/, 2 ,t) : e‘‘=='^®(y, 2 ,t) (33) 

where 2 , t) and 2 , t) are the “tissue” stress 

and strain and C‘'®®“®( 3 /, 2 , t) the tissue stiffness matrix. The 
stiffness matrix is determined at the tissue scale, the explana¬ 
tion is presented in Section 2.1. 

Based on the Bernoulli hypothesis, the strain distribution 
appears to be a plane and remains plane even after deforma¬ 
tion. This is why we can decompose the strain by introducing 
three constants: £ 1 , K 3 and K 2 - 

= £i(i) - K3(i)2/ + K 2 (t)z. (34) 

By introducing this relation into Hooke’s law, we obtain: 

a‘r“®(y,2,i) = C‘’r“"(2/,^,f) {er{t)-K3{t)y^K2{t)z) (35) 

Because we assume the shear force to be null, the stress ten¬ 
sor is reduced to one component: a-fff^'^’^{y,z,t). And with 
Bernoulli hypothesis the strain tensor contains only one com¬ 
ponent. Hence, the stiffness matrix can be replaced by the 
component 2 , t). 

Here we can see that if we determine the strain constants, 
we would know the stress distribution. The mechanical load¬ 
ings, the inputs of this model, allow us to determine the strain. 
Indeed the cross section is supposed to be under a normal force: 
N and a bending moment M here divided in two bending 
moments: My and M^, such as M rh = My y + M^ z. By 
definition of the stress we have the relations: 

N = j a^^--'^%y,z,t)dA (36) 

My = j z-a^^t'^%y,z,t)dA (37) 

M. = -lyalT'^-iy,z,t)dA (38) 

By introducing the static moments of first and second order: 
EA, ESy, ESz, Elyy, EIzz, Elyz, the equations become the 
following constitutive relation: 
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where: 


“ =y 

^ClT'^-{y,z,t)dA 


ESy = J 

^ClT'^%y,z,t)-ydA 


Elyy = j 

''clT'^-{y,z,t)-y^dA 


Es. =y 

^ClT'^-{y,z,t)-zdA 


II 

^ 2. t) ■ yzdA = Eiyz 



^ClT'^-{y,z,t)-z^dA 


If we chose the origin of the coordinate system at the 
normal center (NC) of the cross section, the coupling terms 
between extension and bending vanish since they become null 
by definition of the NC: 

ESy = j 

ClT'^^y,z,t)-ydA = 0 

(40) 


ClT'^'^{y,z,t)-zdA = 0 

(41) 


The constitutive relation can be simplified as: 

’a] \eA 0 0 1 fsi' 

My = 0 Elyy Elyz K2 (42) 

_Mz_ _ 0 Elyz EIzz_ _K3_ 


By definition ESy and ESz with respect to the y - z coordinate 
system are zero. From which the unknown position of the NC 
with respect to the known position of the y — z coordinate 
system can be found: 


Vnc 


ESy 

EA 
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Z'NC 
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To conclude, here is the step-by-step methodology we are 
using to find the stress and strain distribution in the cross 
section: 


1. Localise the normal center (NC) by computing the inte¬ 
grations: EA, ESy, ESz- 

2. Compute the integrations: Elyy, EIzz and Elyz- 

3. Determine the cross-sectional forces: N, My and Mz- 

4. Calculate the cross-sectional deformations: si, K 2 and kq 
based on Eqn. (42). 

5. Find the strain distribution based on the kinematic rela¬ 
tion, Eqn. (34). Here it is important to remember to use 
the coordinate centred in NC. 

6. Find the stress distribution based on Hookes’ law, Eqn. 
(33). 

The initial cross section is extracted from a microradio¬ 
graph, as it is explained in Section 2.4; and the mechanical 
loading is not symmetrical. Hence the position of the NC is 
changing. This is why we need to localise it after each step. 
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Determination of the Normal Force Center NC 

The special location of the origin of the coordinate system for 
which the coupling terms {ESy and ESz) between extension 
and bending become zero is by definition called the normal 
force center NC. The result of this definition is that an axial 
force N which acts at the NC only causes straining and no 
bending. The coupling terms are also referred to as weighted 
static moments or weighted first order moments. To find the 
position of the coordinate system for which the coupling terms 
become zero requires a tool. 

Assume a temporary coordinate system: y — z from which 
the porosity distribution is known. The shift in origin of this 
coordinate system with respect to the y - z coordinate system 
through the unknown NC is denoted with yNC cmd ^ivc- The 
temporary coordinate system can be expressed in terms of 
the y - z coordinate system as: 

y = y + VNC z = z + ZNc 

Hence: 

ESy = I {y,z,t)- ydA (43) 

= I ClT'^^y,z,t)-ydA + yMC f ClT'^^y, z,t}dA 
= ESy + EA ■ ytqc 

ESz = yC*'|=“®(j/,^,t) • zdA (44) 

= f ClT'^^y,z,t)-zdA + z^C I ClT'^^y,z,t)dA 

= ESz + EA ■ zmc 
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